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Abstract. We exhibit a characteristic structure of the class of all regular graphs of degree d 
that stems from the spectra of their adjacency matrices. The structure has a fractal threadlike 
appearance. Points with coordinates given by the mean and variance of the exponentials of 
graph eigenvalues cluster around a line segment that we call a filar. Zooming-in reveals that 
this cluster splits into smaller segments (filars) labeled by the number of triangles in graphs. 
Further zooming-in shows that the smaller filars split into subfilars labelled by the number of 
quadrangles in graphs, etc. We call this fractal structure, discovered in a numerical experiment, 
a multifilar structure. We also provide a mathematical explanation of this phenomenon based 
on the Ihara-Selberg trace formula, and compute the coordinates and slopes of all filars in terms 
of Bessel functions of the first kind. 
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1. A numerical experiment 

For the sake of simplicity we will pay our attention mainly to cubic graphs (or, in other 
words, to regular graphs of degree d = 3). This assumption is not restrictive since all 
our considerations remain valid for regular graphs of degree d > 3 with obvious minor 
modifications. Moreover, in a certain sense cubic graphs are the generic regular graphs 
(see e.g. Greenlaw and Petreschi [3jjE 

So let us consider the set of all regular cubic graphs with n vertices. They can be 
conveniently enumerated using the GENREG program of Markus Meringer [6] . Each graph 
is completely determined by its adjacency matrix, which is symmetric. Its spectrum (the 
set of eigenvalues) is real and lies on the segment [—3, 3]. For each graph it can be found 
numerically. In the interests of statistical analysis, we might want to take the means and 
variances of each set of eigenvalues. However, since the diagonal entries of the adjacency 
matrices are zero (graphs contain no loops), the eigenvalues sum to zero. In order to 
produce results with some variation, and originally motivated by solving systems of linear 
first order differential equations, we take the exponential of the eigenvalues before finding 
their mean and variance. As a final modification, this time motivated by the authors' 



1 Actually, cubic graphs are generic in a wider sense: any graph can be made cubic by a small pertur- 
bation that blows up vertices into small circles. 
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Fig. 1. Plots of mean versus variance for the exponential of the eigenvalues of the doubly- 
stochastic matrices associated with all regular cubic graphs with various numbers of vertices. 



interest in Markov processes, we replace the adjacency matrix A by the related doubly 
stochastic matrix |A The theory of Markov chains then states that the probability of 
being at the j'th vertex after a walk of length i in the graph with each edge equally likely 
to be chosen is the jth element of the vector (|^4)* x, where the /cth element of the vector 
x is the probability of starting at the kth. vertex. 
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Fig. 2. Successively zooming in on the n = 14 plot 



Summarizing, we apply the following procedure. For a fixed even n find the adjacency 
matrices of all regular cubic graphs on n vertices. In each case, divide the adjacency ma- 
trix by three, find its eigenvalues, take their exponential, and then find their mean and 
variance. Each cubic graph is then represented by a single dot on a plot of mean versus vari- 
ance. Figure 1 shows the results of applying this procedure with n = 8, 10, 12, 14, 16, 18, 
where the number of regular cubic graphs in each case is 5, 19, 85, 509, 4060, 41301 respec- 
tively. There appears to be a very definite structure in these plots. In each case the data 
appear in distinct clusters that at this scale look like straight line segments with roughly 
the same slope and distance separating them. (In the next section we will derive explicit 
formulas for these slopes and distances.) Due to their form, we would like to name these 
clusters by "filars", whose dictionary meaning is "threadlike object s"q 

An even greater level of structure exists within each filar. Figure 2(a) repeats the results 
for n = 16, and Figure 2(b) zooms in on the leftmost filar. We can see that each filar is 
in fact made up of smaller clusters of approximately straight line segments, all roughly 
parallel and the same distance apart, with a steeper slope than the original one. We shall 
call each of these clusters a subfilar, and Figure 2(c) zooms in on the fourth subfilar from 
the left in Figure 2(b). The structure continues in Figure 2(d), which zooms in on the 



2 This term also recognizes the second author's (JF's) initial investigation of the above phenomenon. 
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Fig. 3. All cubic graphs with ten vertices, as output by GENREG 

5th subsubfilar of Figure 2(c). Since a fractal is defined as a self-similar image, where the 
same structure is evident when magnifying one part of the image, we see that these figures 
obviously enjoy a fractal structure. The larger the number of vertices, the more levels of 
magnification can be undertaken before the number of data points becomes small enough 
for the self-similar structure to be lost. Collectively, we refer to this phenomenon as the 
"multifilar structure" of cubic graphs (or their spectra, to be more precise) . 

Finally, it is worth noting that this behavior is not limited to cubic graphs. Plots for 
quartic graphs (every vertex of degree four) show exactly the same structure. As we will 
see later, the Ihara-Selberg trace formula justifies the presence of such a fractal structure 
for regular graphs of arbitrary degree d. 

Having - to the best of our knowledge - for the first time discovered this property of 
regular graphs, the aim of the following section is to theoretically explain the (multi)filar 
structure. 
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Fig. 4. A reproduction of Figure l(n = 10) with labels replacing data points related to associated 
with graphs in Figure 3 



2. Theoretical justification 



2.1. Ten vertex cubic graphs in detail 



Before we resort to the theory based on the Ihara-Selberg trace formula [HQ], it is in- 
structive to consider the case n = 10 in detail. Figure 3 shows all 19 regular cubic graphs 
with ten vertices, labelled in the order produced by GENREG, and Figure 4 repeats Figure 
1(b), but with labels on the data points indicating graph number. We can see a pattern 
in terms of which graphs are in each filar. Graphs 19,18,16,17,15,14 in the first (leftmost) 
filar have no subcycles of length 3, which we shall call triangles from now on. Graphs 
13,11,8 in the second filar have exactly one triangle, graphs 10,9,6,7,5 in the third filar 
have exactly two triangles, graphs 12,3 in the fourth filar have exactly three triangles, and 
graphs 4,2,1 in the fifth filar have exactly four triangles. Focusing on the first filar, graphs 
19,18,16,17,15,14 have exactly 0,2,3,5,5,6 subcycles of length 4 respectively. Special em- 
phasis should be given to the data points for graphs 17 and 15. They are extremely close 
together, and have the same number of subcycles of lengths 3 and 4 (0 and 5 respectively). 
They only vary in number of subcycles of length 5, numbering and 2. Similarly, graphs 
4 and 2 both have four 3-subcycles, two 4-subcycles, and only vary in the number of 5- 
subcycles (two and four). These observations suggest that membership in a filar structure 
is related to the number of subcycles of various lengths in the original graph - although 
this does not explain why filars approximate straight lines. 
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2.2. Explicit formulas for the mean and variance 

To obtain qualitative and quantitave justification of the phenomenon in question we bring 
in a very explicit version of the Ihara-Selberg trace formula that is due to P. Mnev, cf. 
[TJ, Formula (30). A general form of the Ihara-Selberg trace formula, as well as precise 
definitions can be found in Appendix. For any regular graph G of degree d = q + 1 on n 
vertices we have 



1 V^- g + 1 /^r* d, I 1 VV ^ T (2 Tat) m 

i=l ^ ' 7 fc=l 

Here {Ai, . . . , A n } is the spectrum of the adjacency matrix of £?, 7 runs over the set of 
all (oriented) primitive closed geodesic^ in G, £(7) is the length of 7, and I m (z) is the 
standard notation for the Bessel function of the first kind: 

r=0 K ' 

All lengths £(7) are integers greater than or equal to 3. Let us denote by the number 
of non-oriented primitive closed geodesies of length I in the graph G. The numbers me are 
called the multiplicities of the length spectrum of the graph G, that is, the set of lengths 
of non-oriented primitive closed geodesies in G (the set {777.3, m i, ■ ■ ■} describes the length 
spectrum of the graph in a unique and convenient way). The multiplicities are uniquely 
determined by the eigenvalues of G and are given by explicit formulas (cf. Formula (34) 
in [71 In. For instance, we have 



1 n I { n \ 

-J^Xl m = - f^Af-n(g+l)(2g + l) , 

1=1 \ 8=1 / 



" 3 (, 

etc. 

Now we rewrite Formula (CD) in terms of the multiplicities (since we consider in detail 
only the case of cubic graphs, we also put q = 2): 



n 00 

_ J2 e Ul = J(t) + - V £m,F,(t), (2) 
77 z — ' 77 z — ' 

i=i e=3 
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and 



3 In the context of graphs, a closed geodesic is an oriented closed path of minimal length in its free 
homotopy class. A closed geodesic is called primitive if it is not a multiple of a shorter geodesic. Closed 
geodesies of length 3, 4 and 5 are 3-, 4- and 5-cycles in graphs respectively, whereas closed geodesies of 
length greater than 5 may have self-intersections. See Appendix for details. 

4 Clearly, the length spectrum also determines the eigenvalue spectrum uniquely. 
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k=l 

Note that the factor of 2 at the sum appears because we forget about the orientation of 
geodesies and have to count each one of them twice. The latter series converges very fast 
because of the following well-known asymptotic behavior of I m {z): 

/.« * (5) (3) 



as m — ► oo and < z <C vm+T; see e.g. [9]. 

The closed form expressions for the mean /i and the variance a can now be easily 
extracted from (jSJ). Precisely, we have 



1 n „ oo 

^ = " E eAl/3 = J (V3) + - V £m^(l/3), (4) 
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o (5) 
= J(2/3) + - X) tn t F t (2/3) - /i 2 . 

n £=3 

Substituting (BJ into the last formula and neglecting quadratic terms in Fe that are 
small in view of (j3J), we get 

2 " 



a w ( J(2/3) - J(l/3) 2 ) + - X ^ (*>(2/3) - 2 J(l/3)F/(l/3)) . (6) 

£=1 

Now we are set for explicitly describing the positions of fllars. 



2.3. Coordinates and slopes of filars 

To start with, let us note that the function Fe(t) is positive for t > and decreases very 
rapidly when I grows and t remains fixed. It is easy to check that ^(2/3) — 27(1/3)^(1/3) 
is also positive for any positive integer i > 3 and decreases very fast in i. Therefore, for any 
n all the points with coordinates (//, cr), corresponding to cubic graphs on n vertices, lie 
above and to the right of the initial point (J(l/3), J(2/3) - J(l/3) 2 ) « (1.17455, 0.4217) 
in the mean vs. variance plane. Moreover, we see that t = 3 gives the leading terms in 
both sums in and (jSJ). This means that the points (//, a) accumulate just above the 
line parametrically described by equations 

x= J(l/3) + tF 3 (l/3) « 1.17455 + 0.00653t, 

y= (J(2/3) - J(l/3) 2 ) + i (F 3 (2/3) - 2J(l/3)F 3 (l/3)) « 0.4217 + 0.0462t, 

It should be mentioned that the plots on Figures 1 and 2 are build using the unbiased variance 
sf l _ 1 = ;^TyX)r=i ( e "^ — AH • I n order to make formulas simpler, we consider here the variance = 

i X)"=i ( elsL — ) ■ The difference is insignificant, especially for large enough n. 
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(a) (b) 





Fig. 5. A string of diamonds (a) and a clasp (b) 

where x is the mean and y is the variance coordinates respectively. Note that the slope of 
this line is approximately 7.079. 

Another byproduct of the above considerations is a necessary and sufficient condition 
for two graphs and to belong to the same filar: this happens if and only if the 
multiplicities and np are equal, or, equivalently, G^ and G^ have equal number of 
triangles. The lines these filars approximate are given by parametric equations 

x= J(l/3) + 6n 3 F 3 (l/3)/n + tF 4 (l/3), 

y= (J(2/3) - J(l/3) 2 ) + 6n 3 (F 3 (2/3) - 2 J(l/3)F 3 (l/3)) /n 
+ t(F 4 (2/3)-2J(l/3)F 4 (l/3)). 

The horizontal distance between two filars that contain the points corresponding to G^ 
and G^ is proportional to np — np and is approximately equal to 

W(M%\ F n /o ^3(2/3) -2J(l/3)F 3 (l/3) \ r£-rg 
F 3 (l/3) - ^(l/3) F4(2/3) _ 2J(1/3)F4(1/3) J . 

For n = 12, 14, 16, 18 the approximate horizontal distances between the neighboring filars 
are 0.00181, 0.00155, 0.00136, 0.00121 respectively. As one can easily see on Figure 1, 
filars actually get closer to each other as n gets larger in proportion with 1/n. However, 
the slope of filars is independent of n and is equal to F 4 (2/3)/F 4 (l/3) — 2J(l/3) ~ 15.89. 
All the above agree perfectly with the numerical data plotted in Figures 1 and 2. 

Each filar splits into subfilars labelled by the number n 4 of quadrangles in the corre- 
sponding graphs. These subfilars approximate line segments of slope i 7 5(2/3)/i 7 s(l/3) — 
2J(l/3) ~ 33.36. The horizontal distance between subfilars is measured by increments of 

8 ( F M n\ F M /o ^4(2/3)-2J(l/3)F 4 (l/3) \ 
n ^ 4(1/3) " F5(1/3) F 5 (2/3)-2J(l/3)F 5 (l/3)J " 

One can pursue this kind of analysis for subfilars of any level (or depth). 

Finally, a comment is in order. The first one concerns the number of filars, or, equiva- 
lently, how many distinct values the multiplicity m 3 (= the number of triangles) can attain 
for a regular graph on n vertices. We sketch an argument that gives an upper bound for 
m 3 in terms of n. It can be shown that the maximal number of triangles is achieved in 
planar graphs. By the Euler characteristic formula we have v — e + f — 2, where v, e, / 
are the numbers of vertices, edges and faces of a planar graph. Denote by fk the number 
of its /c-gonal faces. Then, for a <i-regular graph on n vertices v = n, e = nd/2, f = fk- 
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Substituting these expressions into the Euler characteristic formula, we get 

^ ' A:>3 

Say, for d = 3 a careful analysis of this formula shows that m 3 = / 3 < 2[n/4] with the 
exception of K±, the complete graph on 4 vertices. This upper bound is sharp. For n 
divisible by 4 it is achieved by looping a string of diamonds shown in Figure |5](a) (cf . [3] , 
[5]). When n = 2 mod 4, we need to attach a clasp on either end of a string of diamonds, 
as shown in Figure 0(b). In fact, for cubic graphs on n > 8 vertices the multiplicity m 3 
can be any number between and 2[n/4]; the corresponding examples can also be easily 
constructed (e.g. m 3 = for bipartite graphs). 

Appendix. The Ihara-Selberg trace formula 

The famous Selberg trace formula relates the eigenvalue spectrum of the Laplace op- 
erator on a hyperbolic surface to its length spectrum - the collection of lengths of closed 
geodesies counted with multiplicities. An immediate consequence of the Selberg trace for- 
mula is that the eigenvalue spectrum and the length spectrum uniquely determine one 
another. A similar result is valid for regular graphs [I], pQ. To formulate it precisely we 
need to introduce some terminology. 

We consider oriented closed paths in graphs up to cyclic permutations of vertices. An 
elementary homotopy is a transformation of a closed path of the form 

(V i, ...,Vj,...,V n , Vi) f-> (v h Vj, v', Vj..., v n , Vx), 

where Vj and v' are adjacent vertices. Two closed paths are called (freely) homotopic if 
one can be transformed into another by a sequence of elementary homotopies or their 
inverses. The unique shortest representative in a (free) homotopy class of closed paths is 
called a closed geodesic. The length £(7) of a closed geodesic 7 is the number of edges it 
contains; 7 is called primitive if it is not a power of a shorter geodesic. 

Now let G be a regular graph of degree d = q + 1 on n vertices. Denote by A its 
adjacency matrix, and let {Ai > A 2 > . . . A n } be the spectrum of A. Note that Ai = q + 1 
and I Ai| < q + 1. The following result can be found in [TJ: 

Theorem 1. Let h : Z — > C be a sequence of complex numbers such that h(n) = h(—n) 
for all n G Z and 

00 

^|Mn)|g"/ 2 <oo. 

n=l 

Put h{z) = Xl^L-oo h{ n ) z ~ n , the discrete Fourier transform of h(n). Then 

i=l v II 7 fc=l 

where Zi is related to A» by the equation A» = y/q(Zi + z^ 1 ), and 7 runs over the set of all 
primitive closed geodesies in G. 
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Formula §B) follows from ([7j) if we take h(n) = I n (2^Jqi). Then we have h(z) = 
e t^/q(z+z )_ rj^g integrals that enter ([T]) and ([7j) are related to each other by the change 
of variable s = y/q(z + z^ 1 ). 
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